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ABSTRACT 

We extend electromagnetic finite elements based on a variational principle that uses the electro- 
magnetic four-potential as primary variable. The variational principle is extended to include the 
ability to predict a non-linear current distribution within a conductor. The extension of this the- 
ory is first done on a normal conductor and tested on two different problems. In both problems, 
the geometry remains the same, but the material properties are different. The geometry is that of 
a one-dimensional infinite wire. The first problem is merely a linear “control” case used to validate 
the new theory. The second more interesting problem is made up of linear conductors with vary- 
ing conductivities. Both problems perform exceedingly well and predict current densities that are 
accurate to within a few ten-thousandths of a percent of the exact values. The fourth potential is 
then removed, leaving only the magnetic vector potential, and the variational principle is further 
extended to predict magnetic potentials, magnetic fields, the number of charge carriers and the 
current densities within a superconductor. The new element generated by this formulation is then 
tested on a one-dimensional infinite superconducting wire. The element produces good results for 
the mean magnetic field, the vector potential and the number of superconducting charge carriers 
despite a relatively high system condition number. The element did not perform well in predicting 
the current density. Numerical problems inherent to this formulation are explored and possible 
remedies to produce better current predicting finite elements are presented. 
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into complex problems can be made more productive. It centers on the observation that 
some aspects of the problem axe either better understood or less physically relevant than 
others. These aspects may be then temporarily left {done while efforts are concentrated 
on the less developed and/or more physically important aspects. The staged treatment is 
better suited to this approach. 

1.2 Mechanical Elements 

Mechanical elements for this research have been derived using general variational principles 
that decouple the element boundary from the interior thus providing efficient ways to work 
out coupling with non-mechanical fields. The point of departure was previous research into 
the free-formulation variational principles presented in Ref. [3]. A more general formulation 
for the mechanical elements, which includes the assumed natural strain formulation, was 
established and presented in Refs. [4-7]. New representations of thermal fields have not 
been addressed as standard formulations are considered adequate for the coupled-field 
phases of this research. 


2. ELECTROMAGNETIC ELEMENTS 

The development of electromagnetic (EM) finite elements has not received to date the 
same degree of attention given to mechanical and thermal elements. Part of the reason 
is the widespread use of analytical and semianalytical methods in electrical engineering. 
These methods have been highly refined for specialized but important problems such as 
circuits and waveguides. Thus the advantages of finite elements in terms of generality have 
not been enough to counterweight established techniques. Much of the EM finite element 
work to date has been done in England and is well described in the surveys by Davies [S] 
and Trowbridge [9]. The general impression conveyed by these surveys is one of an un- 
settled subject, reminiscent of the early period (1960-1970) of finite elements in structural 
mechanics. A great number of formulations that combine flux, intensity, and scalar po- 
tentials are described with the recommended choice varying according to the application, 
medium involved (polarizable, dielectric, semiconductors, etc.) number of space dimen- 
sions, time-dependent characteristics (static, quasi-static, harmonic or transient) els well 
as other factors of lesser importance. The possibility of a general variational formulation 
has apparently not been recognized. 

In the present work, the derivation of electromagnetic (EM) elements is based on a vari- 
ational formulation that uses the four-potential as primary variable. The electric field 
is represented by a scalar potential and the magnetic field by a vector potential. The 
formulation of this variational principle proceeds along lines previously developed for the 
acoustic fluid problem [10,11]. 

The main advantages of using potentials as primary variables els opposed to the more 
conventional EM finite elements based on intensity and/or flux fields are, in order or 
importance: 

1. Interface discontinuities are automatically taken care of without any special interven- 
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tion. 


2. No approximations are invoked a priori since the general Maxwell equations are used. 

3. The number of degrees of freedom per finite element node is kept modest as the 
problem dimensionality increases. 

4. Coupling with the mechanical and thermal fields, which involves derived fields, can 
be naturally evaluated at the Gauss points at which derivatives of the potentials are 
evaluated. 

The problems which are presented in this paper only have one geometry, that of an infinitely 

long cylindrical wire (see Fig. 1). The problems that are examined are: 

1. Infinite conductor with the same conductivity between elements. 

2. Infinite conductor with different conductivities between elements. 

3. Infinite superconductor. 


3. LINEAR FINITE ELEMENT WITH CURRENT PREDICTION 

In the previous development of electromagnetic finite elements, the current distribution, J, 
has always been assumed to be given. Unfortunately, J for a superconductor is not linear, 
and with the exception of a few special cases, unknown. Since it is not known a priori 
what the current distribution for a superconducting element will be, these elements must 
have the ability to predict the current density J, as well as the magnetic potential A. The 
starting point in developing a superconducting element is to build a linear element with 
current prediction capabilites. Because current density varies from element to element 
in a superconductor, our new linear element should be able to model a current density 
that changes from element to element. This occurs for a linear media when the material 
conductivity, ct, changes from element to element. 

3.1 Finite Element for a Linear Conductor 


The potential energy functional for a linear conductor is given by 

p= /^{£( VxA M e (- v *) 2 - jrA } w 

The constitutive relation (Ohm’s law) is: 

J — —oV<f> (2) 

where o is the material conductivity. 
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Substitution of the (2) into (1) yields 


0 =J v iV {r^ xA ^-rO 2 - 3TA } (3) 

Variation of the above gives 

SU = J dv{^(V x £A) T (V x A) - C (“) T (“) -^J T A- J t £a} (4) 

For a one-dimensional axisymmetric problem where the only non-zero components of A 
and J ar.e A z and this reduces to 

CU = f v dv{±(^)[^)-j;(SJ 1 ) T {J,)-6J,A,-J,6A 1 } (5) 

Integration by parts produces: 




( 6 ) 


This is appropriate if J and a are continuous. If a is discontinuous, i.e., varies from element 
to element, then we need to examine what happens at element boundaries. Since the energy 
functional has variational index m = 0 for J, we may approximate an element current 
density, Jj, as a step function. Maxwell’s equations require that V x E = a -1 V x J = 0. 
The one- dimensional discretized equation is 


' J « J: u(ri 

Q 

- r<-i)| ri = flr i+i^* i+ i^: u ( r *+i “ r «)U 

(7) 

*7' Win 

- rj.i)| r< = °r+i J t i+ A r <+i ~ r »')lr, 

(8) 


_ — 1 7 C — — 1 Jt 

°i J xi — a i+l J z i+ i 

(9) 


where u(r, -r.) and S(rj -r*) are the Heavyside step and Dirac delta functions respectively. 
This boundary condition equation, weighted by a Lagrangian multiplier, is added to the 
original functional. These multipliers will be denoted by Af f , where i is the number of the 
boundary between two adjacent conducting elements. 

The last requirement on this system of equations comes from the law of charge conservation. 
With I being the total current flowing through a surface T, and n the unit normal, this 
law is 
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or in discretized form 


I = 


L 


J • hdT 


( 10 ) 


nvmeJ - 

i- E / J*n*<fr* = o 

• / r* 


(ii) 


Multiplying this equation by the global Lagrangian multiplier A y and adding to the mod- 
ified energy funcional gives, in discretized form 


numel 


T(U riiCI /* - T® \ T / T c % 

C '= E X.^{^( V xN„Ay r (VxN m Ay- f(^) (^)-^N„a:J 

m=l *m ~ 

nimel- 1 n um el - 

+ E A^-’^wr^J+A,^- E / J *- fi * dr 0 ( 12 ) 

1=1 *=1 *' r * 

Next we change / to J 0 -f A//,, where I 0 is the initial current value, Ii is the amount of 
loading current added to the system, and A is the control parameter. The variations of the 
modified energy functional axe then taken with respect to A e , J e , A*., A y , and A to get the 
internal and external force vectors. 


The internal force vector f is 


^ rr numtl m 

m = l m 


1 <9N m T 3N r 


. A e - N r J t ) 


numtl 


numtl— 1 


dr dr 

KUIJtCI * 

- E / <*KUN-A + 

, Jy 

m = 1 m 

nvmcl- 1 

+ E *?>,*’ -<#.)+ E k -1 - »;+i j ;.+.) 

1 = 1 «=1 

nttmc/ - numei - 

-\(-E / <*n) + ('<•- E / 

Si •'!•> Si •'r* 

where vector v includes the degrees of freedom. The external force vector p is 

, X d ' V \T 

p = Aq = ~ x m\ = ~ Xhex - 

where q is the loading vector and e\ t is the vector with all zero components except the 
one associated with the degree of freedom at A y , which is one. 


(13) 


(14) 
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Letting w be the vector of incremental velocities then 


AA Z 
A A z 


w = < 


AJ Zl 

A\h 

AJ Zi 

AXl 2 


i m«l 


AA/ BtlfB<| _j 

AA, 

Since the tangent stiffness matrix, K, is not separable on an element level, KAv is pre- 
sented here rather than the expression for K. The product is 

numtl m - Tqtvj 

KAv= £ Tr*Ai-NSAjt.) 

m=l "Vm ** 

numtl * 

- E / ^(N m AA* m + «T- J Aj; m ) 

— 1 Jv € 

m= 1 m 

numel-1 numcl-1 

+ E A *f. (»r* -<#.)+ E « lAJ i - ^r+W'j 


1=1 


1=1 
nitmel 


numel • mime/ * 

-AA,( E / <“'*)- E / AJ «. <ir * 

Jfc=l * /r * fc=l * /r * 


(15) 


3.2 Nonconducting Element 


For an element outside of' the conductor (for example, free space) J e equals zero, and 
consequently the energy functioned reduces to 


ntime/ * - 

u= E JjV'{±(V*X m Kf(VxN mK )} 

m=l ^ m 


( 16 ) 
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So on an elemental level, in the one- dimensional case, f e , q e , and K e become 


r- / 


dV, 




,dN, 


( 17 ) 


K* 


q* = 0 



(IS) 

(19) 


4. NUMERICAL VALIDATION OF LINEAR ELEMENT 


4.1 The Finite Element Model 


The finite element formulation derived in the previous section has been applied to two test 
problems described below. Both problems Eire treated with one-dimensional axisymmetric 
elements. Each bar element has two end points, one interior node and a common shared 
global node. These nodes are defined by their Eixial position rj. The two end nodes have 
one degree of freedom each corresponding to A z { and A z j. FYom these values the magnetic 
potential components are interpolated with the standard linear shape functions, which 
provide the C° continuity required by the variational formulation. The interior node is 
placed at the center of the element and the global node at the end of the finite element 
mesh. Neither of these nodes carry Einy physical signifigance and are used solely to provide 
the extra degrees of freedom assigned to the two LagrangiEin multipliers and the degree of 
freedom assigned to J*. J e and A* are carried on the center node and \ g is CELrried on the 
common node. Consequently, each element has 2x 1+2+1 = 5 degrees of freedom. 

For the calculation of the element stiffnesses and force vectors, it is assumed that the 
permeability \i and current densities are uniform over the element. The desired stiffness 
matrix and force vector are calculated by numerical quadrature using a 2 point Gauss 
formula. 


4.2 Applying Boundary- Conditions 


The finite element mesh is necessarily terminated at a finite size. For the two test problems, 
the outer radial end of the mesh is defined as the truncation radius r = Rt. The outer 
radial end of the conductor’s mesh is defined as the wire radius r = R w i re . Since current 
is only carried in the conductor, the degrees of freedom for J e between R w i re and Rt are 
constrained to zero. Constraining A to zero at Rt causes both boundary terms in equation 
(6) to disappear. Notice that constraining A to be zero at r = 0 does not remove the outer 
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boundary term in equation (6). This is shown in section 5.1, where A(r=0) is required to 
be zero. 

4.3 Assembly, Solution and Field Recovery 


The components of the element stiffness matrix corresponding to the summations from 1 
to numel in equation 15 axe first calculated and inserted into the master stiffness matrix. 
The components for the remaining two summations are next determined and added to 
the master tangent stiffness matrix. This is done in an element by element fashion until 
the complete master tangent stiffness matrix is assembled. The loading force vector is 
assembled in an element by element fashion following standard finite element technique. 
The boundary conditions are set as explained in the previous section. The modified master 
equations modified for B.C. are processed by a standard symmetric skyline solver, which 
provides the value of the magnetic potential at the mesh nodes, and the mean current 
density over each element. 

The physical quantity of main interest is not the potential, but the magnetic field B$. 
This is obtained by discretizing A as follows. Since A is only a function of r, B# = 
— (3N/3r)A*. This represents the mean value of B over the volume of the element. In 
previous papers [15, 16], extrapolation of this value to the endpoints was tried with little 
success and was not tried here. Consequently, the value obtained for B is plotted as a step 
function over the elements in our figures. 

The ability of the potential formulation to accurately model the discontinuity in the B 
field at a conductor/free space interface has already been established in previous works 
[15, 16]. For this reason, in both test problems fi is set equal to 1 inside the conductor 
and in the free space surrounding it. The first test problem sets all the Ci s to one, and 
the second problem sets the for the element equal to the element number. 


4.4 Problem 1: Equal Conductivities 


The first test problem is identical to that reported in Schuler and Felippa [15] with a 
one-dimensional axisymmetric discretization. As shown in Figure 1, it consists of a wire 
conductor of radius R w i rt transporting a total current I = 1 in the z direction. For all 
of the problems reported here, the elements were assumed to have a unit thickness in the 
z direction. Since none of the desired quantities vary in the z direction this is a valid 
assumption. The radial direction is discretized with N w i rt elements inside the wire and 
Nfne elements outside the wire in free space. The mesh is truncated at a “truncation 
radius” Rt, where the potential A, is set equal to zero. Other boundary conditions are 
set as previously defined. 

The results obtained with Rt = 2 R w i re , N wire = 20, Nf rt t = 20 for the potentials matched 
those generated by our previous linear electro-magnetic finite elements [15, 16]. Because 
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a potential may vary by a constant, it is not necessary for the two curves to lie one upon 
another, only that they have the same shape. Figure 3 shows the curves obtained for the 
exact and computed solutions for the potential and verifies that these curves are almost 
an exact match. Figure 2 shows the analytical and computed solutions for the current 
densities. The result obtained for the computed current density is lower than the true 
value by less them one ten thousandth of a percent, thus providing a check on the element 
calculations. Because these results are so close to the exact solution, they are plotted as 
a series of points, rather than a line, so that they may be distinguised from the exact 
solution. 


4.5 Problem 2: Different Conductivities 


This problem’s geometry, Rt, Rwire , N/ rt t and N w i re are the same as that reported above. 
The only difference is that now the element conductivity is set to the element number. The 
values obtained for the current densities were just as good as those obtained in the first test 
example. Similarly, the values obtained for the potentials and magnetic fields gave results 
with an order of accuracy of previous finite elements. The computed potential varied from 
the analytical potential by an almost constant value, as Figure 6 shows, and the curve for 
the true value of the magnetic field intersected the middle of the top of the “step” of the 
computed solution. The magnetic field was represented as a step function here, unlike our 
results in previous papers. This was done to more accurately portray that the computed 
B field is the mean value for the true B field over the element. The computed value does 
not go to zero at the element boundary (r = 0), as expected, since it is the mean over 
the first element. Close inspection of Figures 4 and 7 reveals that the true value does not 
intersect the middle of the “step” there, showing that as we get closer to the center of the 
conductor, the computed value is higher than it should be. Numerical experiments reveal 
that if a finer and finer mesh is used, that the computed value converges closer to the true 
value, as expected, therefore verifying the validity of this model to accurately calculate A, 
B, and J. 


4.6 Conclusions with Respect to a Current Predicting Element 


The results obtained in the previous two problems show that it is possible to extend 
previously derived finite element formulations for static and quasi-static magnetic fields 
to cases where the current distribution in the element is unknown. This is particularly 
encouraging since this means that it should be possible to solve problems where material 
and geometric nonlinearities preclude a linear current distribution. It also means that 
whereas before a knowledge of how the current was distributed in a conductor was needed, 
with the new formulation, all that is needed is the total current I through the conductor, 
its material properties /i and a, and the conductor geometry. 
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Table 1 Theory Nomenclature 


Symbol 

Quantities 

a, p 

Temperature dependent material parameters 

+ 

Analagous to a wave/position 
function in particle mechanics 

M 2 

Number of superconducting charge carriers 
per unit volume 

r 

Complex conjugate of tp 

q m 

Effective charge on charge carriers 

m* 

Effective mass on charge carriers 

n 

Planck’s constant divided by 2 ir 

A 

Magnetic potential vector 

B 

Total magnetic field 

J 

Current distribution 

A l 

Lagrangian multiplier whose domain 
is over the volume of interest 

F t 

Helmholtz free energy of superconducting state 

Fn 

Helmholtz free energy of normal state 

AF 

F t - F n 


Having successfully extended our formulation to include the ability to predict current den- 
sities, we next attempt to solve the problem of a nonlinear conductor, the superconductor. 


5. SUPERCONDUCTING F.E. WITH CURRENT PREDICTION 

This section presents the basic theory for superconductivity and how the constraints and 
finite elements for one-dimensional superconducting problems were developed. Table 1 
presents the nomenclature for the quantities used in this work. 

5.1 The Helmholtz Free Energy for a Superconductor 


In the general vicinity of the transition or critical temperature for a type I or II supercon- 
ductor, the change in the Helmholtz free energy ( F ) can be approximated as 


A F = F,-F„ = J dv{-oH , + J^| 4 + 5 ij|(-«W-j , A)^| , + iB.H} (20) 
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in S.I. units [17], where the quantities a, (3 and ip are defined in Table 1. The first two terms 
represent a typical Landau expansion of the Helmholtz free energy for a second order phase 
transition. The third term represents the total momentum of the charge carrier. The —ih'V 
term is analogous to the dynamic (kinetic) momentum of a quantum wave- like particle*; 
the q* A term represents the field momentum [13, p. 633; 14, pp. 105-108]. 

Using the identities, B = ^ 0 H, and B = V x A, the last term, which represents the field 
energy, can be replaced by 

ji-(Vx A.) 2 (21) 

The material’s magnetic permeability fi , is set to fi 0 , the value of the permeability in free 
space. This value is chosen because the field energy term represents the magnetic energy 
in free space, and the other three terms are corrections to that energy resulting from 
material and dynamic effects. Unlike a linear material, corrections to the field energy in 
a superconductor cannot be accounted for by an appropriate choice of a constant fi. The 
magnetic permeability within the conductor is now a function of the spatial coordinates 
and is no longer a constant. 

Expanding A F in terms of ip and ip* gives 

A F= f dvi-onPxP* + ±P(rPiP*) 2 + -^-(-ihV T xP-q*A T xP) 

Jy v 2/72* 

(ihVip* - q*Aip*) + ;dj-(V x A) T (V x A) | (22) 

Talcing ip = ipR + ipj, ip* — ipR — ipj, where ipR represents the real part of the order 
parameter and ipj the imaginary part, gives 

AF = J^dV {-a(ip R 2 +ipi 2 ) + %P(rpR 2 + ipl 2 ) 2 

+ + ~ <fA T (ip R + iipi))(ihV(ip R - iipi) - q*A(ip R - iipi)) 

+ ^(Vx A ) T (VxA)} 

= J v dv { -“(V>R 2 + <h 2 ) + + ’pi 1 ) 2 + i (^(V’VrVV’r + V t ^,V4,,) 

+ 2^ A t (VV>rV>/ - VMr) + A a t a(V>r j + V/ 2 )) 

+ i-(V X Af(VxA)} 

(23) 

* A good example is a one-dimensional particle in an infinitely deep energy well. The — ihV 
term in our functional is similar, in theory, to the momentum of the particle in the well. 
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The first variation of A F with respect to ipR is 


«AfWlO = J <*v{ty R (-2a* « + 2/?(V>s 2 + <h^R - + £1a t A^ r ) 


+ 


v ^£v^*+0 aI V)} 

= J dv{Srl> R (~2ai> R + 2/?(V'r 2 + */*)*« - 0V</,,A T + £ r A 7 'Afe+ 

’(^» + pA^,))} + 

(24) 


The variation with respect to t/>/ is 


6AF(64>i) = f dvh^j(- 2a^/ + 2/?(tf>* 2 + V’/ 2 )^/ + ^V^A r + ^-A r AV'/+ 

v(£v^,- 0 a ^*))}+ 

/ r ™ r {(£ v *'-£ A *-M 

(25) 


The boundary terms in the last two equations represent the net momentum exchange 
between the B field and \ip\ at the boundaries. At the inner boundary, we expect the value 
of rf> to approach a constant for a bulk superconductor. Therefore we require that A go to 
zero. At the outer boundary, we require that there be no superconducting flux into free 
space; to ensure this we set tp equal to zero there. 

The variation with respect to A (again in one dimension) is 
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2 

6&F{6A.) = J dv{6A(££(V T <P R <fo - V t M r ) + ^ r A T (V>« 2 + v; 2 )) + 

i(Vx(A) T (VxA)} 

Po J 

= J dv{6A(££(V T tf> R iJ>i - V T M R ) + ^A T (rP R 2 + rpj 7 )+ (26) 

J-(VXVXA^))}+ 

— f dddzS A{r^-\ 

Ho J r i or ) 

Note that A t = { A r , A$, A z }, which reduces to A r = {0, 0, A z } for the one-dimensional 
case and also that dr = dOdz. Constraining the value of A = 0 at r = R w irt causes both 
boundary terms in the previous variation to vanish. At r = R w irci the value of B will be 
the same as the value for B given by a linear conductor with the same current flowing 
through it. This is a direct consequence of Ampere's circuital law, f t B • ds = pi. From 
V x A = B, the second boundary term becomes 


Po 


J dddzSA {r„„> e B J 


Using Ampere' s law, this becomes 


(27) 


L 


deiz6A{—) 


(28) 


This will give a reaction force at r = R w j re (lifter integrating over dT) of -I. The required 
boundary condition is that A r (r=0) is constrained to zero, but doing this will give a 
reaction force at r = 0 of To maintain the same reaction force, a current I is added to 
the boundary at r = 0 and a current -I is added at r = R w i rc • This gives the residual (for 

_4v(r=0) = 0): 


[ dvi — (V x V x A) + - Vt})]tj) R ) + ^-;A(il> R + V>/ 2 )} = 0 (29) 

J y 'Ho m TTl J 

From Maxwell’s laws, it is known that p~* V x V x A = J; substituting this formula into 
the above equation gives the constraint on J, i.e., 

dV { J + ^(V<M>/ - VV>/0h) + ^A(V>* 2 4- rPi 2 )} = 0 (30) 
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For one dimension, multiplying the above equation by the Lagrange multiplier A i and 
adding it, along with the current conservation constraint A,, to the Helmholtz free en- 
ergy give the desired functional for a superconducting finite element. This topic and the 
mathematics will be taken up in a later section. 

5.2 Evaluation of Material Parameters a and ft * 


Deep inside a superconductor, due to screening effects (the Meissner effect), there are no 
fields or gradients. The functional AF’s last terms drop out and the resulting equation is 


AF = -aM 2 + (31) 

Near the second order phase transition, at the critical temperature, the minimum value 
for the free energy occurs when 

= — 2aM + = 0 (32) 


W 2 = Wool 2 = J (33) 

where |V>oo| 2 is the value for the number density of superconducting charge carriers deep 
within the conductor. Substituting |t/>o ©| 2 back into the preceding equation for A F, gives 


A F = 


a 2 a 2 

T + ^8 



(34) 


When the critical field B c is applied, A F = -B 2 /2p 0 . Because of this condition, deep 
inside a superconductor, where no gradients are present, the following approximation to 
A F can be made 


A F = 



B| = 

2 /? ^ fi 0 ~ P 


The work, W, done in setting up a current distribution J [12] is 


(35) 


W = -\ [ 3 T AdV (36) 

Jv 

From London theory [14, p. 84], with A e // equal to the effective London penetration depth, 
the following equation relating J and A is derived 


1 The following has been abstracted from Tinkham’s textbook [14], pp. 105-109. 
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J = - 




Substitution of this expression for J into the equation for W gives the result 


= . ; 2 — / A T A dV 


From Ginzberg-Landau theory [14, p. 107] , the expression for the work done in setting up 
a current density J is defined as being 


w= W aTa ^ 2dv 


If gradients of the order parameter are zero and there are no external fields present, both 
of the above equations are good approximations to W. Equating these two expressions for 
W gives 




_*2 

hr. KM* 


Algebraic manipulation produces 


IV'coI 2 = — 

i i rt* 


/M* 2A e// fi 


Solving for (3 gives 


From before 


a Po<r K) 
P = a ^ 

m* 


B* a 2 


Po P 

Substitution for (3 and more algebra produces the following values for a and /? 
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Equating IV’ooj 2 with n*, which is the number of superconducting electron pairs, we see 
that to be consistent with London theory 

q* = — 2e = twice the electron charge 

m* = 2m = twice the electron mass 


5.3 Analysis of the London Type Superconductor 


The solution of the discretized superconductor follows the equilibrium curve for A F. This 
curve has at least one critical point, at I = 0. To make numerical experiments proceed, a 
good approximation to the solution vector is needed to move the solution off this point. The 
London type superconductor is such an approximation and is used in this work to advance 
the solution. For a London type superconductor, the gradient terms are assumed to be 
small, and ^>r 2 +^i 2 goes to \ipoo\ 2 over the whole conductor. Making these approximations 
to the previously derived Ginzburg-Landau equations produces 

J = — —“AIV’ool 2 (45) 

m 

Performing some algebra and vector mechanics gives us the differential equations for B 

*2 

Vx J = -2— ’ Vx A|V>oo| 2 ( 46 ) 

m* 


, tf * 2 

— V x (V x (V x A)) = - 2 — V x A|^oo| 2 (47) 

Uo m 


V x (V x B) 
In one dimension, this becomes 


m 


x l/f 


B 


d 2 B e I dBg 

dr 2 ^ r dr 



Letting B 6 j\\ ff = U$ (r), u 2 = 1 and r/A e // = x produces 


&Ue(x) I dUejx) 

dx 1 x dx 



Uo (x) = 0 


(48) 


(49) 


(50) 


This is just the modified Bessel equation whose solution is U$(x) = dl u {x) + b/C»(x). 
J„ and K v are the modified Bessel functions, and a and b are constants dependent upon 
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the boundary conditions. For the solution to remain bounded as it approaches 0, b = 0, 
because fC„ ( x ) — ► oo as x — *■ 0. Using this information, the solution becomes 


U, (x) = all (*) = V, = all 


At r = R w irc, Be = Xy. 0 Ii/2T:R w i rt . Substitution gives 


(51) 

(52) 


a 


h toh 


2nR v 


• X lff 





(53) 


Using the identity V x A = B, the relations A x (r/A e //) = —aXljjI 0 (r/X e ff), and J z 
= aX e ff/n 0 lo ( r/X e ff ) can be derived. The previous formula for J predicts the nodal or 
positional value of J. The mean value for J x is desired rather than this formula because 
any finite element developed in this work predicts the mean of J over an element. Using r,- 
and rj to represent the inner and outermost nodal positions of an element, J z is integrated 
over the volume of the element and divided by this volume to produce the mean value: 


f_ (*Q xdx ^ e // ^ (j0 


J xdx 


Vo 


" =2a ^kfLM 




r j Me// 


• r fM«// 


(54) 
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5.4 One Dimensional Superconducting Finite Element 


For the one-dimensional axisymmetric case, the discretized energy functional, minus the 
boundary terms of equation (26), is expressed as 


nume/ * 

V- £ / iJCl-atofNjN.to. + faMN.k) 

n=l 

+ iW/jNjN.*/. + 

[h 2 (ihlv NjVN„^ /n + 0*rVNjVN„^ s „) 

+ V«ArNl(^*IVNlN.*. - NjN„^ s „) 

+ j* , ArNrN.A.(«,rNj , N.^,. + ^ inin.**.)] 


+ 


+ A e , 


Lm 
2 


+ ijATNr^/ITNrN.^. + ^InTn.^j,.) + J 


+ J_(v x N„A n ) r (V x N.A.) 

Zfl 0 


(55) 


nume/ 

+ Ap |/ 0 + a/l — 

n=l 

Variation of the above with respect to rpR produces the following residuals on an elemental 
level 



r*„. = j <JV'*^ R eT {-oN , 'WJ i + 2 ^N t (^; t N t N^J + <i>! i 7 'N , 'N,^)NV’R+ 

— [VVN t VN</>J ! + q'hlv N r N^5 - NVN^;)NA‘+ 

m* L 

q * 2 N r N^(A e 7 N r NA*)] + 

\\ [££ (VN t NV»/ - N r VN tpj) + 2^-N r N^ A e l ) (56) 

1 m* m J J 
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For ip i, -we get 

r *,• = f dV e 6rP 1 eT {-o i N T NrP t I + 2/?N r (^) T N T N0} + ^ T N T N^)NV>)+ 

-L [n 2 VN r VNV’/ - 9*ft(VN T N^ - NVNV>h)NA c + 

g* 2 N r NV»/(A eT N r NA*)] + 

A£ [-^ (VN r N^ - N t VN^) + 2^-N r NV>/ A e ] } 

For A we get 

r A « = [ dV e 6A eT \ n t ( tp e R T VN T 'Nxpj — xpj T 'VN T Nrp R )^—+ 

Jye. { TTi 

^-N T NA £ (V’/ r N r NV>j + V’r T N 7 ’N^k)+ 

^VISFVNA* + A e L N T (^(V-J r N r N^ + V>R T N r N^))} 

For variation with respect to the J £, s, we get 

rj« = f dV'SJ'Xi - [ dT e 63 e X g 

Jv • Jr* 

Variation of the A^’s gives 

r A « = J dV6\* L { ££ (ij> e R T VN T Nrl> e j - +fV N r N^) + 


^(^ T N T Nt/>R + # r N r Ntff)NA e + J £ + 

+ A/ l 


/**&+££} 


And the variation of \ g (for an element) produces 


SX 


/ Jo + A h _ r 
\ TlUTTlcl Jr *« / 


(57) 


(58) 


(59) 


(60) 


(61) 
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Making the substitutions 


= N il>} 

n =n^ 

A e = NA e 

gives for the internal force vector 


d* e , 

Or 

dr 

ML 

dr 


= VrJ>} 

= VA e 


where 


( 

f V’J 


r = 


f A - 

f J‘ 


f *i = / i/< ^'{N 7 '(2»^(- a + ^(tf + n 2 ))-0^M' + -'i)+ 

+ 2*1)) + + 0¥}(a(' + *1))} 


U; = J v /V’{N T (2^’ l (-c, + 0(<lf + n i )) + ^^(A’ + y L )+ 

+ 2Ai)) + VN r (£^ _ + Ai))} 


i 

m 


(62) 


(63) 


(64) 


r A* 


= dV<{ N^(^(^^-^n)+^(^+n 2 )M e +Ai))+VN T ^ 0 - 1 } 

(65) 


fj« = [ dV'Xi - [ <£T e A, 
Jv j r« 


( 66 ) 


f ^ = X. dy, {S^-^n)+^m s +^M‘+J'+ / r ^4> (67) 


fA, = — - 2 — : - [ dT e y 

TlUTTLCl jp* 


( 68 ) 
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The loading vector q is 


-oM 


dr 

dX 



( 69 ) 


m 


Taking the second variation with respect to the v/s produces the tangent stiffness matrix 


3SW = Jv iV ‘ {A ' T {^ N M^ r N r N« + ^ r N r NVi) + ^VN r VN}«A' 

(70) 

dIW R = L dVWT {" T W T V T ™ ~ N t N)^+ 

2 £!n i '(NA' + Ai)(N«)N}^J (71) 

dJJW, = J v . <tV'SA‘T{nT(rn T VN T N - AWVN)£? + 

2^N r (NA' + Ai)(N05)N}^5 (72) 


a 2 (7 

dA e d\ e L 



dV € 6 A eT 



N r (v>) T N T Nv-j + ^ t n t n^) 


(73) 


a 2 i7 

5A e aj e “ 0 


(74) 


d 2 u 

dA e dX g ~ 0 


(75) 


m 
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( 76 ) 


d 2 U 

d'lpftdipft 



dV e 6xl> e R T { N r N(-2a + 2/?(V>/ 7 N r NV’| + 3^H T N r N^)) + 
2 

^•VN r VN + N r N (A * r N T ) (NA e + 2A |) 


d 2 E 7 

dxl> R dip} 


J dV e 6iP% T {N T N4p(xj>} T N T NxP%)+^ (NA e +A£) (V T N T N-N r VN) }<5^ 

(77) 


a 2 t/ 

drp R dX e L 


J <fy e ^ T {0(VN r NV>)-N T VNV>|) +N r 2^N^NA e }a^ (78) 


d 2 U 

d 2 £T 

dip e R d\ g 


(79) 


(80) 


0 2 t/ 

dipjdxfrj 


J <fy e <Stf>) r {N r N(-2a + 2 / 3 (V»h T N t N ^ + 3V’j T N r N^)) + 


^ 2 T 
— VN r VN + 

m* 



N T N(A er N T )(NA e + 2X1 



(81) 


d*U 

dip}dX e L 


J <iy e (5^ r |0(N r VN^- VN r N^) +N T 2^-NV>/NA e }(5A e L (82) 


d 2 U 

dipjdJ* 

&U 

dxpJdXg 


&U 

dX e L dX e L 


= 0 


d 2 U 

dX e L dJ e 



dV e 6\ e L 6J* 


(83) 


(84) 


(85) 


( 86 ) 
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( 87 ) 


d 2 u 

d\ldX g 

d 2 u 

dJ e dJ e 


( 88 ) 


d 2 U 
dJ'dX g 



dT'SJ'SXg 


(89) 


d 2 U 

dXgdXg 


(90) 


Using the same notation as in the internal force vector f*, the tangent stiffness matrix may 
be represented as 


'Pr ’K-'Pr'P) 

K e w e = 

symm. 

with 


K VJA- 

K *k*l 

®2xl 

k v;A- 

K^ja* 

02x1 

k a*a* 

K A‘a. 

02x1 


0 

x *l* 


0 



K W* = 

dV e |n t 2a + 20 (tyf + 3$k 2 ) + ^—A* (A e + 2X € L )j N + ^ VN r VN J 

(92) 


K Wi = 

j vt dV e |N r (4^5^ e ft )N + ^(A e + X e L ) (VN T N - N T VN)| (93) 

V = 

L dv ‘ ( 4 '5 vn7 ' n - ^Vn) + 2 £ (■*• + nN r N} (94) 
K *l»l = X. iV ‘ {0 (™ T *5 - N r ^) + N r ( 2 ^*M e ) } (93) 
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= 

f dV e j N r (-2ot + 2/5 + 3 (A e + 2A£)^ N + ~VN r VN| 

Jv« l V m m (9 6 ) 

K *,*A* = 

j <JF'/0(^N 7 'N-nVN , 'N^+2^(>t t + Ai)5';N r N} ^ 

**!»! = X. - nVN r ) + N T ^tM*) } («) 


K A . A < = / <fv* | ^-N T N (*; 2 + n 2 ) + X vnTvn } (99) 

K AWv .^'{£ Nr (*' ,2 + **)} (100) 

K x .y=J V 'JV‘ ( 101 ) 

*J-> = / <fT* (102) 

* ./r* 


6. NUMERICAL EXPERIMENTS 

The finite element formulation described in the preceding section has been applied to 
the solution of a one-dimensional axisymmetric superconducting wire. This problem was 
treated with a bar-like element. Each element contains two end nodes, one centroidal node 
and a common global node as in the previously described linear element. These nodes are 
defined by their radial position . Unlike the linear problem, this element contains three 
degrees of freedom on the t\yo end nodes, namely t/vj , V’ij i A j. Standard linear shape 

functions are used to interpolate these quantities across an element and provide the C 
continuity required by the variational formulation. The centroidal node carries no physical 
signifigance, and is used to provide the extra degrees of freedom for the elemental current 
densities and Lagrangian multipliers, and J e . The common global node is assigned the 
degree of freedom for X g . This gives each element a total of2x3+2xl + l = 9 degrees of 
freedom. 

For the calculation of the tangent stiffness, internal force vector and the loading vector, 
the permeability is a constant fi 0 . The London penetration depth and intrinsic Pippard 
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coherence length for aluminum were used to calculate the critical magnetic field B c and 
the effective London penetration depth A,// [13, p.337]. The subroutine that generates 
the elements in the coding, calculates the material parameters a and and assigns them 
to each superconducting element. The tangent stiffness matrix, internal force vector and 
external force vector are calculated by numerical quadrature using a 2 point Gauss formula. 

6.1 Applying Boundary Conditions 


The finite element mesh is terminated at Rj, as in the linear conductor. To ensure that 
no superconducting flux can cross the conductor’s outer edge into free space, the degrees 
of freedom corresponding to J, xpR and xpi between Rwire &nd Rt are set equal to zero. At 
r = R w irei TpR and V>/ a re also set to zero. By setting these values to zero, the boundary 
terms in equations (24) and (25) are required to go to zero. This choice also ensures that 
no superconducting flux can cross either of the boundary surfaces. At r = 0, the value for 
A is also set equal to zero, as discussed in section 5.1. This choice for A requires that |^| 
become constant as the inner boundary of a superconductor is approached. 

6.2 Assembly and Solution 


This section will present the assembly techniques and the solution method used on the test 
problem. Table 2 presents the nomenclature used for quantities of interest. 

The tangent stiffness, internal force and loading vectors are assembled following standard 
finite element techniques. The tangent stiffness matrix K, is stored using a symmetric 
skyline storage scheme, and then modified for boundary conditions. 

Since this is a nonlinear problem, nonlinear techniques must be used to solve for the dis- 
placements of the variables A*, %l> e R , ipj, J e , X e L , and A^. An incremental scheme is employed 
to advance the solution along the equilibrium curve by making a good approximaton to 
the exact solution. An iterative technique is then used to converge to the exact solution. 
The constraint used to limit the distance travelled along the equilibrium curve is arclength 
control, i.e. , |As n | — where A s n is the distance along the curve, and l n is the length 
of the increment, an input variable. The formula implemented in our coding for arclength 
control is listed below in its scaled and unsealed forms. Also listed are the formulas for 
the vectors a and g. The superposed tilde represents a scaled quantity. 

Unsealed Form 


As n = |wjAv n 

Jn 

+ AA„| - /„ 

( 103 ) 

a T = w n // n 

1 

9 ~ In 

( 104 ) 
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Table 2 Solution Nomenclature 


Symbol 

Quantities 

K„ 

Master tangent stiffness matrix at increment n 

q 

Loading vector 

r„ 

Residual vector at increment n 

f„ 

External force vector at increment n 

v„ 

Solution vector at increment n 

Av n 

Vn + 1 V n 

A n 

Control parameter (A=0, zero load, A=l, full load) 

w 

Incremental velocity vector = K“*q 

In 

Input arclength constraint 

|As| 

Distance along equilibrium path 

fn 

y/l + w^w 

s 

Scaling matrix 

C 

Constraint equation 

a 

dc/dv 

9 

dc/dX 

< 

Solution vector at step n and iteration k 

d 

v *+! — 
v n y n 


X at step n and iteration k 

V 

A‘ +1 - A* 

J fact & fact ••• 

first scaling factor for J, A, etc. 

3 scale » Gscalty 

second scaling factor for J, A, etc. 

J/iomopeneouj * ®/iomo0eneoii*i ... 

third scaling factor for J, A, etc. 


Scaled Form 

As n = i- | wjAv n + AA n | - l n (105) 

Jn 

- a T = w „//„ g = 1 //„ (106) 

Av ss SAv, q = S -1 q, K = S _1 KS _1 , 

w = Sw, / = y/l + w T S 2 w = y/l + w T w (107) 

A forward Euler integration scheme was used to predict Av n , A„+i, and v n+1 . Appropriate 
formulas are listed below. 
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(108) 


AA„ = A„+i - A„ = 7 L s!'sn(q T w), 

Jn 

Av„ = K“ 1 q n AA n = w„AA„ 
v„+i = v n + Av n 

A conventional Newton- Raphson technique was implemented for the corrector.This gives 
the linear system 

(.“ ?]{;) 

Since this augmented system is not symmetric, 
are solved for instead 


= -{c} (109) 

the two linear symmetric systems below 


to get 


Kd r = — r Kd, = q 


c + a T d r j 

V = ; — frr-, d = d r + T}d g 


g + aTd, 


which will finally give 


vj +1 = vj + d, A £ +1 = Aj + 7) 

We then iterated until the 2-norm was less than a given tolerance e. 


6.3 Field Recovery 


( 110 ) 


( 111 ) 


( 112 ) 


The primary quantity of interest is not any of the variables in the solution vector, but the 
derived quantity B. This quantity can be computed by using one of two methods. The 
first is to use the identity B = V x A. Discretization in one dimension produces 

B# - -fr A ‘ = (113) 

where Aj and A* are the values for A at the outer and innermost nodes of the element 
respectively, and l e is the element length. 

The second method is to use Ampere’s law. Discretization in one dimension gives 
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2tt r „ ’ ^enclosed — * ^ Jf{ r l+1 ~ r l) (H4) 

where r n is the position of the outermost node of the element on which B is being evaluated. 
The innermost and outermost nodes of elements that are interior to r n are r,- and r, +] , 
respectively. 

The advantage to the second method is that it can be used to predict the B field at the 
end nodes of the elements. The first method merely computes the mean of the field over 
the element. In the test problem, both methods are used to compute the B field and 
compared. 


6.4 Numerical Nightmares in Solution Method 


Starting the Solution Process 

At A = 0, the tangent stiffness matrix is singular if I Q = 0. K has a rank deficiency of 1 
at this critical point. The eigenvector has only one component at the degree of freedom 
corresponding to \ g . Multiplying the eigenvector by q T gives us a non-zero value and tells 
us that this is a limit point. To move off of the limit point, initially a random perturbation 
technique was employed. All of the components of w were perturbed randomly in an effort 
to move the solution off of the limit point. This technique is not recommended because 
there is no guarantee that the approximate solution will be close enough to the equilibrium 
curve to iterate back onto it. Another failing is that the solution may progress back to the 
limit point because the perturbation pushed the guessed solution below the limit point on 
the equilibrium path. For the test problem studied here, this method failed approximately 
eighty percent of the time. 

Instead, to get off the limit point, first all of the Lagrangian multipliers are constrained 
to zero, and then the results are fed into the corrector. When the results are fed into the 
corrector, the multipliers no longer are constrained to zero. The coding then iterates to 
the exact solution. 

If the solution proceeds toward another critical point, the analytical solution for a London 
(extreme type I) superconductor is inserted into v, and is allowed to iterate from that 
point towards the solution. "The London solutions contain the modified Bessel functions Jo 
and Jj. The values for Jo and Ji are calculated using polynomial expansions [18] that are 
accurate to approximately 10 -7 and then inserted into the v vector of the calling program. 

Scaling of Variables 

The arclength constraint is particularly sensitive to inhomogeneous physical dimensions 
in v. The different variables in v in this problem all have different physical dimensions. 
To improve stability of the solution method and to reduce the condition number of the 
master tangent stiffness matrix, a scaling was done in three parts on the variables in v. 
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The first part consists of scaling the variables to have a homogeneous physical dimension. 
The second part consists of scaling the new variables to reduce the condition number on 
the tangent stiffness matrix. The second scaling is done because although the dimensions 
in v may be homogeneous, one variable may only move 10 -6 units per incremental step 
while another may move 10 6 units. The third scaling is a homogeneous scaling of all of 
the variables. This scaling can be used to reduce the order of magnitude of /„ so that the 
Newton- Raphson iterations becomes more stable. The total scaling factor for each variable 
is then represented by 


for example, 


S 


variable 


t;aria&/e /oc 10 vg ™ Me *~ , « 

VCiri able homo geneo us 


(115) 


_ j/.JO'— 1 - 

3 homogeneoui 


(116) 


In the first part of the scaling, K is given units of volume, and then w T is scaled so 
that it has units of energy per unit volume. With L,M,T , and Q representing units 
of distance, mass, time and charge respectively, this means that v should have units of 
M 1 / 2 L 1/2 /T. This also means that A and \\ have units of ML/TQ , xpi and xpR have units 
of 1/X 3 / 2 , J has units of Q/TL 2 , and \ g has units of M 1 / 2 !?/ 2 /T. For v = Sv, a/„ c has 
units of Q/M 1 / 2 L 1 ! 2 , jf ac has units of M 1/,2 L 5 ^ 2 /Q, rpj ac has units of M 1 / 2 L 2 /T, A^y ac 
has units of Q/M 1 ^ 2 L 1 / 2 , and A s ^ oc has units of Q/M l l 2 L 2 l 2 . In previous experiments 

with linear electromagnetic finite elements, it was found that scaling A by \i 0 gave the 
desired units and stability, so aj ac is set to Using the London approximation to the 

Ginzberg-Landau equations gives jf ac = | — fio l ^ 2 \ 2 t jj\. Numerical experiments produced 

q*B c \ e ffr mean /m* i / 2 as a reasonable scaling value for V’/«c» with r meon being the mean 
value of the lengths of all of the elements. The value of A is set to the same value as 
af ac because in their unsealed forms A and X C L have the same physical dimensions. Finally 
A g f ac is set to fio R-w\rt “ a result of numerical experiments. 

The second part of the scaling is done by performing numerical experiments. Two different 
variables can be monitored to determine what the values for a #ca / e , j, ca U i 'Ptcait, A i aca i c > 
and A fliCa/e should be. If v is monitored, the value of A g§calt is changed until v at the 
degree of freedom corresponding to \ g is of the order of magnitude of 10°. Then j sca u, 
A' L acale , dscaie > and upscale are set in this order, by reducing, at the corresponding degree 
of freedom, the largest value of v to an order of magnitude of 10°. 

The second variable that can be monitored is 1/da, where da is the value of the i th element 

rp ' 

in the D matrix from the L DL decomposition of the master tangent stiffness matrix. In 
this scheme, the largest value of d,j, for the degrees of freedom corresponding to xp is 
reduced to an order of magnitude of 10° by adjusting upscale • Then a Jca / e , A £ Jca , e , and 
jscale are adjusted in the same manner. Finally A Jca j e is adjusted in the same manner, 
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but then more numerical experiments are performed by adjusting it up and down one unit. 
While it is adjusted up and down, the condition number of the matrix is monitored. The 
direction which produces the lower condition number is chosen, and A g scale is adjusted 
until the lowest condition number is achieved. For the numerical experiments performed 
here, the condition number was estimated by using a random perturbation technique. 

No matter which variable is chosen, da or v, sometimes these methods fail because K 
becomes singular. They usually fail because of a bad value for A £, aca j e . If this happens 
it is recommended that the same procedures are followed, but get as close to reducing 
values to an order of 10° as possible. The last scaling value that produces a non-singular 
K is the one that should be used. It is recommended that the values for Vacate, A Sjsca;e , 
fljca/e, ^L SC ait' j scale be variables that are interactive input because sometimes many 

numerical experiments must be performed to achieve the best results. It is also recom- 
mended that the coding be set up such that when any of the second scaling factors axe 
incremented one positive unit, the corresponding degrees of freedom in v are scaled down 
by a power of 10, e.g., 10 10 becomes 10° when 10 is the input scaling factor. This helps 
to eliminate confusion as to what the input scaling factor should be. If the variable da is 
monitored, the da's that correspond to the degrees of freedom for ip and A will be scaled 
in a similar manner, but the remaining degrees of freedom will move by a power of 10 2 
when the scaling factor is incremented one positive unit, e.g., 10 10 becomes 10° when 5 is 
input as a scaling factor. Both of these options have been implemented in the coding for 
the test problem presented here. 

The third and final scaling factor is set by monitoring /„. The idea here is to reduce /„ 
to an order of approximately one, as can be seen from equation (108). Unfortunately, 
sometimes doing this can raise the system’s conditon number to an unacceptable value. 
This factor should only be used to reduce a system’s condition number or to stabilize an 
already unstable solution process. 

6.5 Test Problem 


For our test problem, we use a finite element mesh similar to the one described in the section 
on linear conductors. The big difference between this problem and a linear conductor 
problem is that the charge carriers act like a fluid flow in a pipe with resistance. Most of 
the interesting physics occurs in a thin boundary type layer at the conductor/free space 
interface. Because of this phenomena, we generate a regular mesh of Nbuik elements in the 
interior of the conductor, and a geometric mesh of Nboundary elements in the boundary 
layer. The mesh ratio, Nbuik, and N boundary are input values. Because the element for a 
free space magnetic field has been validated many times before, no elements were generated 
external to the conductor. 

The leading term for the London solutions to J, A, and B is exp(xj - x w i re ). Using this 
information, the depth of the boundary layer has been set to 5.x 575. xA e // . This value 
was chosen for the boundary layer depth because it can be determined where the value 
for exp(xj — XtvtVe), where Xj is between Xboundary and x w i re , I s no l ess than 10 -250 . The 


31 



choice of 10 -250 ensures that machine underflows do not occur if the London solutions are 
calculated to move off of a critical point. Setting the boundary depth to 575. x X e ff does 
not allow the London solutions to go beyond machine precision over the whole boundary 
layer, but does increase the condition number on the tangent stiffness matrix. Numerical 
experiments showed that setting the boundary layer to five times this depth reduced the 
tangent stiffness condition number considerably. 

The values of -12, -3, -5 -15, and -13 were chosen for a, ca u , ttaa/e, j, C a/e, »caU, and 
A g lca i e respectively. These values were chosen by monitoring v and adjusting the appro- 
priate scaling factor as described in the previous subsection. The value for the third scaling 
factor was 10 6 and was chosen by monitoring the value for /„ and the tangent stiffness 
condition estimate. 

All graphs presented here are for A = .00103 with a tangent stiffness condition estimate 
of approximately 4 x 10 6 . Figures 8 and 9 display the results for \ip\ 2 normalized by 
\iplc,. Figure 8 shows the behavior for \^normaliztd\ over t ^ ie whole mesh and Figure 9 the 
behavior in the boundary layer. Since aluminum is a Type I superconductor, these results 
match very well with expected physical behavior. The boundary conditions are seen to 
match well in that at both boundaries the slope for \ip n ormalized\ 2 approaches zero (see 
equations 24,25). 

Results for J are displayed in Figures 10 and 11. Figure 10 shows the behavior for J over the 
mesh, and Figure 11 the behavior in the boundary layer. The effects of the high condition 
number become apparent in Figure 11. Although the general behavior as exhibited in 
Figure 10 physically follows that of a Type I superconductor, inside the boundary layer J 
jumps up and down instead of increasing exponentially. The reason for the high condition 
number also becomes apparent. For Type I and II superconductors, B is excluded from the 
better part of the conductor. To match this physical situation, A and J must become zero 
for values of r between the center of the conductor and the boundary layer. The degrees of 
freedom in this region for A and J in the finite element formulation become approximately 
equal to zero and cause the condition number for the tangent stiffness matrix to rise. 

The results for the B field are shown in Figures 12 and 13. The two different methods for 
calculating B a s discussed in section 6.3 are displayed. Again Figure 12 shows behavior 
over the whole mesh, while Figure 13 shows behavior over the boundary layer. Surprisingly, 
in view of the high condition number, the finite difference method yields reasonable values 
for the mean value of B. The values of B calculated by using Ampere’s Circuital Law do 
not perform well because the values for J are not accurate. 

The finite element model also produced results for A, but these results are not presented 
here because the quantities of primary interest in this application are |V’| 2 ) J and B. 
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7. CONCLUSIONS AND FUTURE WORK 


The results obtained with the one-dimensional superconducting element axe encouraging. 
Although the values obtained for J were not very good, the model performed well ehough to 
give us reasonable physical answers for B, A and \rj >\ 2 . They show that the variational for- 
mulation can provide good results for electromagnetic quantities inside a superconductor. 
The maun difficulty is to reduce the tangent stiffness condition number so that reasonable 
values for J may be obtained! A reformulation of the terms involving A and J is suggested, 
since these terms go to zero over the majority of the model. A way to do this is to express 
B as the sum of an equivalent linear magnetic field, and the magnetic field due to the 
material nonlineaxities (the magnetization field /i 0 M). This gives the following equation 
for the last term of equation 22 

— (V X A + B applied) (^7 X A + B applied) (U7) 

where V x A now is equal to /i 0 M. J is similarly split into J applied and Jmatena/- This 
formulation should help to alleviate most of the conditioning problem because for the bulk 
of the conductor, J applied = — J material and Aappiied = ~A ma teriab The solution that will 
be sought over this area, instead of being zero, will be just the negative of the linear 
solution. Any further conditioning problems will be addressed with any of the previous 
scaling techniques, or a diagonal scaling (Jacobi preconditioner) will be applied. 

After the conditioning problem has been solved, the next step in the development of this 
element will be to add thermocoupling effects. Good semi-analytical approximations that 
relate temperature change to B c and A e // are found in Tinkham[ 14]. These two variables, 
and their temperature dependence, directly effect the material parameters a and /?. This 
temperature dependency will be added to the current coding along with an expression for 
the energy change due to temperature variation. If time permits, the final step in this 
research will be to extend the current formulation to two dimensions. The problem that 
will be studied will be a one dimensional infinite superconductor whose critical current 
has been exceeded. The superconducting flux varies in the radial and axial directions. 
Although finite difference formulations of this problem have yielded good results, none of 
them have been able to adequately match experimental data. The suspected problem has 
been that these formulations do not capture effects due to Joule heating of the wire. A 
two-dimensional superconducting element with thermal degrees of freedom appears to be 
the perfect vehicle to test that hypothesis. 
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J vs. radial distance, all elements equal conductivities. Computed 
solution plotted as points to highlight accuracy of method in matching exact solution. 
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Figure 3. 

A vs. radial distance, all elements equal conductivities. 
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Figure 5. 

J vs. radial distance, element conductivities equal to element number. 

Computed solution plotted as points to highlight accuracy of method in matching exact solution. 
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Figure 9. 

|^’IV|V’oo| 2 vs. radial distance, values in boundary layer plotted. 
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